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Abstract. Detailed information on potentially suitable habitats and forecasted alterations thereof under climate change 
scenarios are critical for the conservation planning of endangered taxa, in particular those with small distribution ranges. 
The Huu Lien Tiger Gecko, Goniurosaurus huuliensis, is a is a micro-endemic species in northern Vietnam. The species is 
listed in the IUCN Red List as Critically Endangered and in CITES Appendix II due to habitat loss and overexploitation 
for the international pet trade. Climate change has been globally acknowledged to impact on many species and it likely has 
negative influences on G. huuliensis. In this study, an ensemble modeling technique is employed, trained with climate and 
vegetation cover conditions, to identify the contemporary potential distribution of this species and assess its alterations 
under different climate change scenarios. Our predictions suggest that the current potential distribution of G. huuliensis 
mostly covers the known sites of occurrence and their surroundings. These areas will narrow significantly and/or shift 
towards higher latitudes under novel climate conditions as can be expected according to future IPCC scenarios. To safe- 
guard in-situ populations of G. huuliensis in the context of the potentially severe impacts, we provide a core refugia map 
that identifies key regions for priority conservation measures, including Lang Son Province and small sites in Bac Giang 
and Thai Nguyen Provinces, northern Vietnam. We highly recommend that the Huu Lien Nature Reserve be selected as a 
“centre” to kick-off conservation actions for the target species. 
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Introduction 


The Huu Lien Tiger Gecko, Goniurosaurus huuliensis OR- 
Lov, RYABOv, NGUYEN, NGUYEN & Ho, 2008, which is one 
of five known Goniurosaurus species in Vietnam, has only 
been recorded from the type locality in the Huu Lien Na- 
ture Reserve (NR), Lang Son Province (ORLOV et al. 2008, 
NGUYEN et al. 2009, NGUYEN 2011). This endemic species is 
found in the evergreen forest on isolated karst formations 
at altitudes from 300 to 370 m a.s.l. (ORLOV et al. 2008). 
Although discovered recently, G. huuliensis was assessed 
as Critically Endangered (CR) in the IUCN Red List of 
Threatened Species (NGUYEN 2018), and this species was 


included in CITES (the Convention on International Trade 
in Endangered Species of Wild Fauna and Flora) Appen- 
dix II and the Vietnam Government's Decree No. 06/2019/ 
ND-CP (Group IIB) in 2019 (NGo et al. 2019b). 

Recent investigations provided convincing evidence for 
the endangered status of G. huuliensis. In particular, this 
habitat specialist has, to date, only been recorded from a 
very small karstic range within the Huu Lien Nature Re- 
serve, with the estimated extent of occurrence (EOO) be- 
ing 30 km? (ORLOV et al. 2008, NGUYEN et al. 2009, NGUY- 
EN 2018). Similar to other Tiger Geckos, the effective pop- 
ulation status of G. huuliensis is expected to be extremely 
small (Neo et al. 2016, 2019a, b). Furthermore, wild Tiger 
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Vulnerability of Goniurosaurus huuliensis to climate change 


Geckos, including G. huuliensis, have been locally over- 
harvested for the international pet trade (NGo et al. 2019a, 
b). These severe impacts may be driving this species to the 
brink of extinction. 

Climate change is a global threat to biodiversity and 
ecosystems, which drives large-scale shifts in species dis- 
tribution, and may lead to a decline of species abundance, 
and even extirpation or extinction of many terrestrial or- 
ganisms (HuGGETT 2004, THOMAS et al. 2004, PARMESAN 
2006, Huey et al. 2009, MONASTERSKY 2014, PrmM et al. 
2014, URBAN 2015, TAYLOR & KUMAR 2016, MARKLE & Ko- 
ZAK 2018, WEISKOPF et al. 2020). Besides direct anthropo- 
genic impacts, climate change may also affect G. huuliensis 
in a negative manner. On the one hand, as a poikilothermic 
lizard, basic physiological functions, such as locomotion, 
growth and reproduction are determined mainly by envi- 
ronmental conditions (ARAUjo et al. 2006, FITZGERALD 
et al. 2018, NGo et al. 2019a, VICENTE et al. 2019). On the 
other hand, a species within a small distribution range and 
a habitat specialist being dependent on specific ecological 
niches will be less capable of responding to novel environ- 
mental conditions as a result of climate change (HUGGETT 
2004, ORLOV et al. 2008, SANDEL et al. 2011, Primm et al. 
2014, MARKLE & Kozak 2018). Thus, this target species is 
considered to be particularly susceptible to climate change. 

To date, no in-situ conservation action has yet been 
implemented to protect wild populations of G. huulien- 
sis. General conservation plans have already been pro- 
posed for all Tiger Geckos (NGo et al. 2019b). However, 
in-situ conservation measures have been only conducted 
to safeguard another threatened congener, Goniurosaurus 
catbaensis, from northern Vietnam, after NGo et al. (2016, 
2019a) provided detailed insights into its population sta- 
tus and utilized microhabitat, and Le et al. (2017) used a 
Maxent approach to predict its potential distribution for 
conservation. 

A considerable gap between designing general con- 
servation plans and practical application of conservation 
measures can limit the efficiency level of conservation ef- 
forts (DUDLEY & PARISH 2006). In fact, conservationists 
normally encounter an obstacle for priority conservation 
efforts when trying to identify where potentially suitable 
habitats and core refugia are (BAUMGARTNER et al. 2018), 
even though all populations should be protected in the 
same manner. Recently, the development of species distri- 
bution models (SDMs), based on species’ geographic co- 
ordinates and their environmental data, have been able 
to predict the potential distributions of species (GUISAN 
& ZIMMERMANN 2000, GUISAN & THUILLER 2005). Sev- 
eral techniques to compute SDMs have been recently ap- 
plied for some rare lizards in Vietnam (VAN SCHINGEN et 
al. 2016, LE et al. 2017, Neo et al. 2019a, Neo et al. 2021). 
However, the modelling of endangered or rare species has 
its limitations when only a few recorded occurrences are 
combined with many predictor variables, which may lead 
to the model overfitting (BREINER et al. 2015). Recently, the 
ensembles of small models (ESMs) approach, a novel strat- 
egy of fitting SDMs, has been employed to overcome the 


limitations of using a single SDM technique and improve 
the qualitative outcome of predictions where only a small 
number of presence coordinates is available (BREINER et al. 
2015, COLA et al. 2017). 

The present study aims to predict the potential distri- 
bution of G. huuliensis based on geographic locations and 
environmental variables (including climate and vegetation 
data) under both current conditions and different future 
scenarios. The ESMs technique is employed to compute 
the models and to project the ensemble through space and 
time. We hypothesize that the currently suitable area of 
G. huuliensis will shrink significantly under the impacts of 
climate change. We also intend to identify potential core 
refugia of G. huuliensis, based on the simulated outcomes, 
to improve the efficacy of in-situ conservation measures. 


Materials and methods 
Study area 


The study area (within 20-24°N and 104-109°E) with al- 
titudes ranging from 1 to 2,139 m a.s.l. (Fig. 1) was select- 
ed as the background area, encompassing all known dis- 
tribution ranges of the five Goniurosaurus species occur- 
ring in northern Vietnam (GrisMER et al. 1999, Vu et al. 
2006, ORLOV et al. 2008, ZIEGLER et al. 2008, NGUYEN et 
al. 2009, NGUYEN 2011). Goniurosaurus huuliensis is known 
only from the karst forest at its type locality in the Huu 
Lien NR, Lang Son Province, northern Vietnam (ORLOV et 
al. 2008). Further populations of G. huuliensis are expect- 
ed to be found in surrounding areas, including all districts 
of Lang Son Province and adjoining provinces (e.g., Bac 
Giang, Cao Bang, Quang Ninh and Bac Giang Provinces), 
where we directly observed karst habitats resembling those 
known from the type locality. 


Data collection 


Occurrence data of G. huuliensis were compiled from lit- 
erature, direct observations, and interviews with local peo- 
ple (Ortov et al. 2008). We conducted field surveys in the 
Huu Lien NR in April and August 2019, and from June to 
July 2020. An extensive initiative interviewing local peo- 
ple was carried out in July 2020 in all districts of Lang Son 
and surrounding provinces (Bac Giang, Cao Bang, Quang 
Ninh, Thai Nguyen) in order to possibly detect as yet unre- 
corded populations of G. huuliensis. 

Coordinates of each captured individual were recorded 
with a GPS Garmin 64, but will be shared upon request. A 
total of 80 occurrence records in the WGS84 projection of 
G. huuliensis were initially collected. However, several re- 
cords were removed to reduce the density of recorded lo- 
cations by using spatial filtering in the packages “dismo” 
and “sp” in R v 3.1.2 (R Core Team 2018) due to duplicates 
or near-duplicates for some localities. Only one occurrence 
locality was randomly selected within each 1-km square. 
Spatial filtering can improve the quality of prediction mod- 
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els by decreasing geographical bias, autocorrelation effects, 
and uncertainty (VELOZ 2009, RADOSAVLJEVIC & ANDER- 
SON 2014). A total of 27 representative occurrence points 
were finally used for modelling. 

For environmental data, we selected bioclimatic and 
vegetation cover variables. Nineteen climatic variables for 
current (averages between 1960-1990, version 1.4) and fu- 
ture conditions from the Coupled Model Intercomparison 
Project Phase 5 (CMIPs5) were obtained from Worldclim 
(https://www.worldclim.org/) (Hijmans et al. 2005). To 
predict the future potential distribution of the target spe- 
cies, we selected future climatic data from two general cir- 
culation models (GCM), including Community Climate 
System Model version 4 (CCSM4) (GENT et al. 2011) and 
Beijing Climate Center - Climate System Model 1-1 (BCC- 
CSM 1-1) (Wu et al. 2014) as expected for the 2050s (aver- 
age 2041-2060) and the 2070s (average 2061-2080). Two 
climate change scenarios of representative concentration 
pathways (RCPs): RCP 4.5 and RCP 8.5, representing in- 
termediate and the most severe levels of accumulation of 
greenhouse gas concentrations in the future climate, were 
used for each model, respectively (Moss et al. 2010, VAN 
VUUREN et al. 2011). For vegetation cover data, we extracted 
the Normalized Difference Vegetation Index (NDVI) and 
Enhanced Vegetation Index (EVI) from within the study 
area using Moderate Resolution Imaging Spectroradio- 
meter (MODIS) sensor data from five years, 2015 to 2019 
(https://earthexplorer.usgs.gov/). The mean, maximum, 
minimum, median, range, standard deviations (STD), and 
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mean, maximum, minimum of warmest and coldest quar- 
ters for the enhanced vegetation index (EVI) and normal- 
ized difference vegetation index (NDVI) were generated 
using Quantum GIS (QGIS Version 3.12.0, Development 
Team. 2020, available at http://qgis.osgeo.org [downloaded 
on 25 March 2020]). 

Eleven variables were eventually kept for the predic- 
tion approaches based on the selection previously made 
by Noo et al. (2021) for another tiger gecko, G. lichtenfel- 
deri. Particularly, six climatic variables (i.e., Bio-2: Mean 
Diurnal Temperature Range, Bio-3: Isothermality, Bio-9: 
Mean Temperature of Driest Quarter, Bio-15: Precipita- 
tion Seasonality, Bio-18: Precipitation of Warmest Quar- 
ter, Bio-19: Precipitation of Coldest Quarter) were used to 
train the climatic model, while five remaining variables 
of vegetation cover (e.g., NDVI of Mean Coldest Quarter, 
NDVI of Minimum Coldest Quarter, NDVI of Minimum 
Warmest Quarter, NDVI of STD, and EVI of Range) were 
used for a separate model based on the vegetation struc- 
ture. The climate and vegetation cover variables were used 
separately to model the potential distribution, in order to 
account for both broad-scale climatic factors delimiting 
the species’ potential distribution and fine-scale habitat 
availability within this potential distribution. Future sce- 
narios of potential vegetation structure are not available 
for limiting these variables. This procedure allowed us to 
disentangle the potential impact of climate change from 
those of land-use changes, which are largely unpredict- 
able. 
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Figure 1. Map of study site in northern Vietnam and southern China, including the selected buffer representing a radius of 50 km 
around occurrence points (violet circles) of Goniurosaurus huuliensis. 
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Species distribution modelling 


We fitted the ensembles of small models with six model- 
ling techniques: artificial neural network (ANN), clas- 
sification tree analysis (CTA), generalized linear models 
(GLM), generalized additive models (GAM), generalized 
boosting regression models (GBM), and maximum entro- 
py modelling (Maxent). All models were calibrated with 
the R package “ecospat” ver. 3.1 (BREINER et al. 2015, COLA 
et al. 2017). To build an Ensemble of Small Models, we used 
subsets of the environmental predictors to create bivariate 
predictor combinations. In particular, six climate and five 
vegetation predictors resulted in 15 and 10 bivariate mod- 
els (BiVa) in each of six selected modelling techniques, re- 
spectively (step-1, Supplementary Fig. Sı). In step-2, we 
evaluated each of these bivariate models using a cross-val- 
idated model evaluation index as a Somers’ D (BREINER 
et al. 2015). Bivariate models with a Somers’ D lower than 
o (i.e., AUC < 0.5) were set to zero and not used to build 
ESMs. One ESM prediction was built for each modelling 
technique in step-3 (ie, ANN-ESM based on bivariate 
ANN models, CTA-ESM on bivariate CTA models, etc.). 
In the final steps (Steps 4 and 5), we built a final ensemble 
prediction (EP-ESM) - the fitted ESM by averaging across 
these ESMs, again using Somers’ D as a weighting factor 
(Supplementary Fig. S1) (BREINER et al. 2015). 

All models were calibrated with presence-only data and 
10,000 random background points selected within a buff- 
er area. As the buffer we defined an area with a radius of 
50 km around the occurrence points using the R packag- 
es “dismo” and “raster”. The fitted ESM was also projected 
to areas beyond the selected buffer within the study area 
to show up alternative potential refugia in the context of 
future climate change. Fifteen-fold cross-validation, with 
subsets of 70% training data and 30% test data, was used 
to evaluate the models. We used three adjustment indices 
to evaluate model performance, including the AUC of the 
Receiver Operating Characteristic Curve to discriminate 
presence from absence (or background), which was then 
used to build a Somers’ D (D = 2 x (AUC - 0.5)) (FIELDING 
& BELL 1997, LOBO et al. 2008, ELITH & GRAHAM 2009), 
TSS is particularly (TSS = sensitivity + specificity - 1) use- 
ful for the modelling of rare species and can be used to 
compare different modelling techniques (ALLOUCHE et al. 
2006). In these indices, values closer to 1.0 indicate better 
model performance (BREINER et al. 2015). 

To assess the predictive capabilities of our ensemble 
model projections, multivariate environmental similar- 
ity surface (MESS) analyses were used to quantify poten- 
tial extrapolation errors (ELITH et al. 2010). The respective 
MESS functions of the “ecospat” package were used (CoLA 
et al. 2017). The MESS analysis compares the environmen- 
tal similarity of any given grid cell within the study site to 
a reference set of grid cells of chosen predictor variables. 
It is used to identify extrapolation in areas with novel en- 
vironmental conditions beyond the training range of the 
models, as is indicated by negative values (ELITH et al. 
2010, 2011). In this study, similarity/novelty was classified 


into four levels: “< -10” (high extrapolation), “-10-0” (low 
extrapolation - margin), “0-10” (low interpolation - mar- 
gin) and “> 10” (high interpolation). We also performed 
MESS analyses to create maps for each future climate sce- 
nario and to evaluate the alteration of novel climatic condi- 
tions from the reference conditions under the current cli- 
mate model. 


Core refugia 


The core refugia for G. huuliensis were identified within 
highly suitable areas in terms of climate, which have val- 
ues above the 10% training presence threshold (high oc- 
currence probability). This is a stricter criterion for con- 
verting the probability maps to binary maps with smaller 
suitable habitats. We also identified buffer refugia with the 
least-suitable environmental conditions with values above 
the minimum presence threshold (intermediate occur- 
rence probability). They were all combined with suitable 
areas as identified in the model based on vegetation indi- 
ces. To identify priority areas and assess the effectiveness 
of the protected areas for safeguarding the target species, 
we collected all shapefiles of nature reserves within the oc- 
currence region of G. huuliensis from https://www.protect- 
edplanet.net/ (accessed in June 2019) (Fig. 1). The identi- 
fied environmental refugia for the target species were af- 
terwards layered with the protected area. 


Results 


As a result of our extensive interviews, we documented the 
occurrence of G. huuliensis in four districts of Lang Son 
Province (including the Huu Lien NR) and another district 
in Thai Nguyen Province, northern Vietnam. According to 
local people, all animals were found in the evergreen forest 
on karst formations. They were observed moving or resting 
on rocky cliffs or hiding in crevices. 

Regarding the model evaluation, Maxent.Phillips-ESM 
and ANN-ESM were the best models, showing the high- 
est mean values of adjustment indices with low variation, 
while GLM-ESM exhibited the lowest value of these indices 
with high variation (Supplementary Figs S2—-S3). Compared. 
to these, the fitted ensemble models (EP-ESMs) of climate 
and vegetation improved the accuracy of prediction as in- 
dicated with higher average values of adjustment indices. In 
particular, TSS, AUC and Somers’D in the fitted EP-ESMs 
performed well with average values of 0.80, 0.92 and 0.83 
in the climate model, respectively (Supplementary Fig. S2). 
These values were relatively lower in the vegetation model 
(0.71, 0.88 and 0.75, respectively) (Supplementary Fig. S3). 
All selected climate and vegetation variables significantly 
influenced the prediction and their contribution rates were 
approximately equal (Supplementary Tables S1-Sz2). 

Under current climatic conditions, the potential distri- 
bution of G. huuliensis covers mainly the sites of occur- 
rence and their surroundings within the buffer area. We 
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recorded a few additional potential sites beyond the buffer 
area in Vinh Phuc Province, northern Vietnam, and in bor- 
der areas between China and Vietnam (Cao Bang Province) 
(Fig. 2A). The suitable climate range covers a total area of 
63,773 km’, of which approximately 42.5% (27,124 km?) 
represent highly suitable habitat (Fig. 2B). However, within 
the buffer area, the suitable habitat only covers 9,577 km? 
(about 15% of the total suitable area). With regard to the 
vegetation model, suitable habitats of G. huuliensis are 
relatively scattered within the study site, but they are also 
gathered mainly around the occurrence sites. In particu- 
lar, in the buffer area, the suitable area of vegetation cov- 
ers 6,874 km? within the intermediate suitable site of cli- 
mate, accounting for 71.8% (Figs 2A-2B-4). In terms of the 
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MESS analysis, the buffer area is filled by a large area of 
interpolated habitat, covering 9,620 km? (approximately 
65.0%) in the current climate model and 8,250 km? (ap- 
proximately 56%) in the vegetation model, whereas we only 
recorded a very small area of extrapolation in both these 
models (Figs 2C-2D). 

With respect to future predictions, in general, the cli- 
matically suitable area of each future scenario fluctuated 
depending on RCPs (e.g., 4.5 and 8.5) and periods (by the 
2050s and 2070s). Under future climate scenarios, each 
predicted a significantly smaller area of the potential dis- 
tribution compared to the contemporary one and implied 
a shift towards higher latitudes in the future (Figs 2-3-4). 
However, we recorded a novel case of the CCSM4 - 
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Figure 2. Projected potential distribution for Goniurosaurus huuliensis following (A) the climate model; (B) the vegetation model 
(blue to red colours indicate higher suitability). Multi-Environment Similarity Surface (MESS) map of novel habitat following (C) the 
climate model; (D) the vegetation model (teal colour represents high interpolation habitat, aqua colour - low interpolation, coral 
colour - low extrapolation, brown colour - high extrapolation). 
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RCP 4.5 scenario by the 2050s, covering 10,311 km? of the 
suitable area within the buffer. Afterwards, this suitable 
area contracted considerably under the same conditions 
by the 2070s (Fig. 4). The range contraction is expected to 
be largest under the RCP8.5 scenario by the 2070s. In par- 
ticular, the CCSM4 circulation model suggests that the po- 
tential distribution may cover only 3,507 km? (representing 
24% suitable area of the current model) (Figs 3D-4), while 
in the BCC-CSM 1-1 model, no highly suitable habitat for 
G. huuliensis was documented (o km?) and the interme- 
diate suitable habitat within the buffer was only 2,555 km? 
(Figs 3D-4). Given the MESS projections, interpolated 
habitats within the current climate-based potential distri- 
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bution of G. huuliensis will be gradually replaced by high 
novelty (extrapolated) habitats in the future (Supplemen- 
tary Figs S4—Ss). It is noteworthy that novel climatic con- 
ditions were recorded in most of the entire buffer area un- 
der the RCP 8.5 scenarios by the 2070s (Supplementary 
Figs S4A4-S4B4). 


Discussion 
Prediction models 


Comparing the adjustment indices of TSS, AUC and Som- 
ers’ D, we confirmed that the ensemble SDM significantly 
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Figure 3. Projected potential distribution for Goniurosaurus huuliensis under different future climate scenarios (blue to red colours 
indicate higher suitability). The average map of two circulation models of BCC_CSM-1-1 and CCSM4 in the scenario of (A), RCP-4.5 
by the 2050s; (B) RCP-4.5 by the 2070s; (C) RCP-8.5 by the 2050s; (D) RCP-8.5 by the 2070s. 
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improved the accuracy of model prediction for rare spe- 
cies. In particular, the average values of these indices in the 
fitted ensemble models (EP-ESMs) of climate and vegeta- 
tion were generally higher than for a single modelling tech- 
nique (Supplementary Figs S2-S3). 

All selected climate and vegetation variables influenced 
the predictions, and their contributions were approxi- 
mately equal in each ensemble model (Supplementary Ta- 
bles S1-S2). However, the specific biological relationships 
between the selected environmental variables (macrocli- 
mate and vegetation cover) and the species remain unclear. 
Studies on two Tiger Gecko congeners, G. catbaensis and 
G. lichtenfelderi, recently described their microhabitat use 
and identified vegetation (canopy coverage) and ambient 
parameters (temperature, humidity, weather condition) as 
the most important characteristics (NGo et al. 2019a, NGO 
et al. in press). Therefore, it is undeniable that the selected 
environmental variables, which might play important roles 
in the natural history of G. huuliensis, constrain the current 
distribution range and even have effects on evolutionary 
processes (see also PYRON & BURBRINK 2009, RODDER & 
ENGLER 2011, ZHANG et al. 2014, RATO et al. 2015, HEIDARI 
2019, NOGUEIRA et al. 2019, SHEU et al. 2020). 

In this study, we found in the current MESS map that 
novel climatic conditions will impact on most sites of oc- 
currence of the other three Vietnamese Tiger Geckos 
(namely G. catbaensis, G. lichtenfelderi, and G. luii), and 
that the potential distribution of G. huuliensis does not or 
only slightly overlaps their occurrence sites (ORLOV et al. 
2008, ZIEGLER et al. 2008, NGUYEN et al. 2009, NGUYEN 
2011, NGO et al. 2019a, NGO et al. 2021). The climate-based 
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potential distribution of the granite-adapted G. lichtenfel- 
deri, as simulated by the Maxent approach, showed that it 
does not overlap with the extent of occurrence of any karst- 
dwelling Tiger Gecko in Vietnam, apart from G. catbaensis 
(Noo et al. 2021). We assume that macroclimatic niche di- 
vergence may have a central role in explaining the diversi- 
fication in Goniurosaurus (GRINNELL 1917, RODDER & EN- 
GLER 2011, FISHER-REID et al. 2012). 


Implications for conservation 


Under the future potential impacts of climate change, the 
potential distribution of G. huuliensis tends to narrow and 
shift towards higher latitudes in the border area between 
China and Vietnam. Similar changes were also previously 
predicted by using the Maxent approach for G. catbaensis 
and G. lichtenfelderi (Le et al. 2017, Neo et al. 2021). Since 
all Goniurosaurus species are locally-distributed specialists 
adapted to unique ecological niches (NGUYEN et al. 2009, 
NGUYEN 2011, NGO et al. 2016, 2019a, 2021, NGO et al. in 
press), we only identified currently suitable areas within 
the selected buffer that can serve as refugia for sustainable 
in-situ conservation rather than areas suitable for coloniza- 
tion. Based on the forecasted climate and vegetation maps, 
we have generated a core-refugia map that identifies key 
regions for priority conservation. The green patches in the 
map (Fig. 5) represent highly suitable habitats with regard 
to both climate and vegetation, mostly in areas in Lang Son 
Province and at small sites in Bac Giang and Thai Nguyen 
Provinces, northern Vietnam. The yellow areas indicating 
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Figure 4. Projected habitable areas of suitable categories for Goniurosaurus huuliensis under different climate conditions of current 
and future scenarios (blue line represents the total areas of intermediate suitability within the study site; orange line the total areas 
of high suitability within the study site; grey line the areas of intermediate suitability in the buffer; yellow line the areas of high suit- 


ability in the buffer). 
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Figure 5. Proposed refugia throughout the range of Goniurosaurus huuliensis under the current model (green areas indicate core 
habitats of suitable vegetation within yellow areas of high climatic suitability and red lines cover all climatically suitable areas in the 


buffer area). 


the intermediate level of climate suitability should be con- 
sidered buffers of the core refugia (Fig. 5). 

There are three nature reserves (namely Bac Me, Huu 
Lien and Than Sa - Phuong Hoang) within the selected 
green patches (Fig. 5). However, the target species has been 
found only in the Huu Lien NR. A large area of evergreen 
forest in the Huu Lien NR is rigorously protected by local 
people and rangers, and wild populations of G. huulien- 
sis may potentially be found in previously overlooked but 
suitable areas there (NGo et al. pers. obs). Thus, the Huu 
Lien NR plays a very important role and should be con- 
sidered a “centre” for the initiation of in-situ conservation 
programs for G. huuliensis. 

However, our interviews with local rangers showed that 
dealers are usually contacted by local hunters from the 
Huu Lien NR with offers of desired quantities and pric- 
es for each wild specimen. Captured animals were then 
handed over to the dealers in places outside the nature re- 
serve to evade stringent inspections by local rangers. As a 
consequence of these illegal actions, large numbers of wild 
animals were collected, leading to a potentially significant 
decline of the wild population of G. huuliensis (Nco et al. 
2019a, b). Furthermore, we observed other anthropogenic 
impacts on the karst habitat of G. huuliensis, such as ex- 
panding road construction and illegal timber logging. 

Thus, we recommend the establishment of a species 
and habitat conservation area for the threatened spe- 
cies G. huuliensis. The core-protection area should be se- 
lected from within the green patches in the core-refugia 
map with the “Huu Lien NR-Center” (Fig. 5). However, 
this suggestion may not be the best option to protect the 


species effectively. In fact, our selection of core refugia is 
only based on the outcome of our climatic and vegetation 
models. To enhance the practical applicability of in-situ 
conservation measures, conservationists should be pro- 
vided with detailed biological information on the target 
species, such as population status, habitat requirements, 
and threats as well. At the moment, to mitigate the an- 
thropogenic impacts on wild populations and their hab- 
itat, we strongly propose that monitoring of the illegal 
trade and protecting natural forests be intensified, as well 
as community education on the value of biodiversity be 
enhanced. 

Two other lizards, Gekko canhi and Scincella apraefron- 
talis, were recently discovered in the Huu Lien NR (Ncuy- 
EN et al. 2010, RÖSLER et al. 2010, NGo et al. pers. obs). We 
assume that these endemic species, being sympatric with 
G. huuliensis, are likewise being negatively affected in the 
context of climate change. Therefore, these lizard species 
will also benefit from improved conservation measures in 
the region. 
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Supplementary Figure S1. Schematic of the Ensemble Small 
Models (ESM) for each selected modelling technique. 


Supplementary Figure S2. Performance of seven Ensemble of 
Small Models of climate according to adjustment indices of TSS, 
AUC and Somers D. 


Supplementary Figure S3. Performance of seven Ensemble of 
Small Models of vegetation according to adjustment indices of 
TSS, AUC and Somers’D 


Supplementary Figure S4. Multi-Environment Similarity Surface 
(MESS) map of the novel habitat following future circulation 
models. 


Supplementary Figure S5. Predicted areas of novel habitats as per 
MESS analyses under different conditions of current and future 
scenarios. 


Supplementary Table S1. Relative contributions (percentages) of 
climatic variables for ESMs. 


Supplementary Table S2. Relative contributions (percentages) of 
vegetation variables for ESMs. 
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